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Abstract. In this paper we deal with a semilinear hyperbolic chemotaxis model in one space 
dimension evolving on a network, with suitable transmission conditions at nodes. This framework is 
motivated by tissue-engineering scaffolds used for improving wound healing. We introduce a numerical 
scheme, which guarantees global mass densities conservation. Moreover our scheme is able to yield 
a correct approximation of the effects of the source term at equilibrium. Several numerical tests 
are presented to show the behavior of solutions and to discuss the stability and the accuracy of our 
approximation. 
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1. Introduction 

The movement of bacteria, cells or other microorganisms under the effect of a chemical stimulus, represented 
by a chemoattractant, has been widely studied in mathematics in the last two decades, see 21,23,26 , and 
numerous models involving partial differential equations have been proposed. The basic unknowns in these 
chemotactic models are the density of individuals and the concentrations of some chemical attractants. One of 
. — ■ the most considered models is the Patlak-Keller-Segel system [T5], where the evolution of the density of cells is 
' described by a parabolic equation, and the concentration of a chemoattractant is generally given by a parabolic 
5— j ■ or elliptic equation, depending on the different regimes to be described and on authors' choices. The behavior of 
this system is quite well known now: in the one-dimensional case, the solution is always global in time, while in 
two and more dimensions the solutions exist globally in time or blow up according to the size of the initial data. 
However, a drawback of this model is that the diffusion leads to a fast dissipation or an explosive behavior, and 
prevents us to observe intermediate organized structures, like aggregation patterns. 

By contrast, models based on hyperbolic/kinetic equations for the evolution of the density of individuals, 
are characterized by a finite speed of propagation and have registered a growing consideration in the last few 
years [5HTl [T51l2"o] . In such models, the population is divided in compartments depending on the velocity of 
propagation of individuals, giving raise to kinetic type equations, either with continuous or discrete velocities. 
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Here we consider an hyperbolic-parabolic system which arises as a simple model for chemotaxis: 

u t + v x = 0, 

v t + X 2 u x = 4> x u - v, (1.1) 
4> t ~D <j) xx = au- b(f>. 

Such kind of models were originally considered in [27] . and later reconsidered in |llj . They are based on an 
adaptation to the chemotactic case of the so-called hyperbolic heat or Cattaneo or telegraph equation, adding 
a source term accounting for the chemotactic motion in the equation for the flux. The function u is the density 
of cells in the considered medium, v is their averaged flux and cf> denotes the density of chemoattractant. The 
individuals move at a constant speed A, changing their direction along the axis during the time. The positive 
constant D is the diffusion coefficient of the chemoattractant; the positive coefficients a and 6, are respectively 
its production and degradation rates. 

These equations are expected to behave asymptotically as the corresponding parabolic equations, but 
displaying a different and richer transitory regime, and this is what is known to happen at least without 
the chemotactic term. Analytically, these models have been studied in [TB1[T7] and more recently in [T3] , where 
the analytical features were almost completely worked out, at least around constant equilibrium states, where it 
is proved that, at least for the Cauchy problem, the solutions of the hyperbolic and parabolic models are close 
for large times. 

The novelty of this paper is to consider this one dimensional model on a network. More precisely, we consider 
system in the form (|l.ip on each arc of the network, and so we have to consider one set of solutions (u,v,(f) 
for each arc. Functions on different arcs are coupled using suitable transmission conditions on each node of 
the network. Conservation laws or wave equations on networks have already been studied, for example in [8] 
for traffic flows or in 4, 29] for flexible strings distributed along a planar graph. However, here we consider 
different types of transmission conditions, which impose the continuity of the fluxes rather than the continuity 
of the densities. Therefore, in this article, a particular care will be given to the proper setting and the numerical 
approximation of the transmission conditions at nodes, both for the hyperbolic and the parabolic parts of (|1.1[) . 
In particular, some conditions have to be imposed on the approximation of the boundary conditions, in order 
to ensure the conservation of the total mass of the system. Let us also mention that a first analytical study of 
system (jl.ip on a network, coupled through transmission conditions of this type, is carried out in [T2"] . 

The study of this system is motivated by the tissue-engineering research concerning the movement of 
fibroblasts on artificial scaffolds |5ni[2H[2B]j during the process of dermal wound healing. The natural process 
of healing of a damaged tissue occurs through a first phase in which fibroblasts, the stem cells to be in charge 
of to the reparation of dermal tissue, create a new extracellular matrix, essentially made by collagen, and, 
driven by chemotaxis, migrate to fill the wound. In recent years, tissue-engineering research has developed 
some new techniques, which aim at accelerating the wound healing. Actually, cellular migration on an injured 
body zone is stimulated and improved by using artificial scaffolds constituted by a network of crossed polymeric 
threads inserted within the wound, which mimic the extracellular matrix. The fibroblasts's reparation action is 
accelerated, since they already have a support and also they are constrained to move along the network, with 
less degrees of freedom, and it is believed that this approach could be effective in minimizing scarring 28 . 
Therefore, our simple model of chemotaxis on a network, which can be obtained by reducing the kinetic model 
of cell movement on a 3D extracellular matrix proposed in [5] to the case of a network, see [T2], is a good 
candidate for reproducing this configuration: the arcs of the network stand for the fibers of the scaffold and 
the transport equations give the evolution of the density of fibroblasts on each fiber. However, in this paper, 
we only address the numerical aspects of this problem. More direct applications of this framework to the real 
biomedical problem will be explored in future research. As reported in [20], the present understanding of the 
critical biochemical and biophysical parameters that affect cell motility in three-dimensional environments is 
quite limited. Nevertheless, it has been observed that junction interactions affect local directional persistence 
as well as cell speed at and away from the junctions, so providing a new mechanism to control cell motility by 
using the extracellular microstructure. Therefore, mathematical modeling and simulations could play a crucial 
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role in providing a better understanding of these phenomena, and an optimization tool for designing improved 
scaffolds. 

The main focus of this paper is on the construction of an effective numerical scheme for computing the 
solutions to this problem, which is not an easy task, even for the case of a single arc. In that case, non 
constant highly concentrated stationary solutions are expected and schemes which are able to capture these 
large gradients in an accurate way are needed. The main problem is to balance correctly the source term with 
the differential part, in order to avoid an incorrect approximation of the density flux at equilibrium, as first 
observed in [14) . Asymptotic High Order schemes (AHO) were introduced in [25], inspired by |Tj, to deal with 
this kind of inaccuracies. These schemes are based on standard finite differences methods, modified by a suitable 
treatment of the source terms, and they take into account for the behavior of the solutions near non constant 
stationary states. An alternative approach, inspired by the well-balanced methods, has been proposed in [51110] . 
with similar results. However the methods in [25] seem easier to be generalized to the present framework. 

Regarding the problem considered in this paper, the main difficulty is in the discretization of the transmission 
conditions at node, also enforcing global mass conservation at the discrete level. Therefore, in Section [5] we 
explain some analytical properties of problem (|1.1[) . with a particular emphasis on boundary and transmission 
conditions. Section [3] is devoted to the numerical approximation of the problem based on a AHO scheme with 
a suitable discretization of the transmission and boundary conditions ensuring the mass conservation. In the 
present paper, we have chosen to consider only the second order version of the scheme, which is enough for 
our purposes, but it is easy to adapt also the third order schemes proposed in [25] . Remark that here, unlike 
the single interval case, we are forced, for any given time step, to fix the space step on each arc using relation 
(|3.14[) introduced in Section 3, to obtain consistency on the boundary. Numerical tests (not shown) confirm the 
necessity of this supplementary constraint. 

Finally, in Section [4] we report some numerical experiments, to show the behavior and the stability of our 
scheme. A special attention is given to the stability of the scheme near nodes and the correct behavior of the 
approximation for large times and near asymptotic states. It has to be mentioned that during this research 
we observed, in contrast with what happens for the diffusive models, the appearance of blow-up phenomena 
even for data of relative moderated size. Even if, up to now, there are no rigorous results, which can help to 
decide if these singular events are really occurring, or they are just a numerical artifact, our close investigation 
in Subsection 4.3 gives a strong indication towards the first alternative. 



Let us define a network or a connected graph G = (AT, A), as composed of two finite sets, a set of P nodes 
(or vertices) M and a set of N arcs (or edges) A, such that an arc connects a pair of nodes. Since arcs are 
bidirectional the graph is non-oriented, but we need to fix an artificial orientation in order to fix a sign to the 
velocities. The network is therefore composed of "oriented" arcs and there are two different types of intervals 
at a node p G N : incoming ones - the set of these intervals is denoted by I p - and outgoing ones - whose set 
is denoted by O p . For example, on the network depicted in Figure [T] 1,26/ and 3, 4 e O. We will also denote 
in the following by I out and O ou t the set of the arcs incoming or outgoing from the outer boundaries. The ./V 
arcs of the network are parametrized as intervals a; = [0, Li], i — X,...,N, and for an incoming arc, Li is the 
abscissa of the node, whereas it is for an outgoing arc. 

2.1. Evolution equations for the problem 

We consider system (ll.ip on each arc and rewrite it in diagonal variables for its hyperbolic part by setting 



Here u + and u are the Riemann invariants of the system and u + (resp. u ) denotes the density of cells following 
the orientation of the arc (resp. the density of cells going in the opposite direction). This transformation is 



2. Analytical background 




(2.1) 
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Figure 1. First example of network. 



inverted by u = u + + u and v = X(u + — u ), and yields: 



u+ + Xut = — ({<f> x - X)u+ + {<j> x + X)u ) , 

Ut - Xu x = {{<t>x - + (4>x + \)U~) , 

— n f ti~ >r 



t - D(f) xx = a(u + + u ) - b<j). 



(2.2) 



We can also denote by = (4> x =F A) the turning rates (namely the probabilities of cells to change direction) 
and a(u + + u~) — b<f) represents the production and degradation of the chemoattractant. We assume that all 
the cells are moving along an arc with the same velocity A (in modulus), which may depend however on the 
characteristics of the arc. For the moment, we omitted the indexes related to the arc number since no confusion 
was possible. From now on, however, we need to distinguish the quantities on different arcs and we denote by 
uf,Ui , Vi and fa the values of the corresponding variables on the i-th arc. On the outer boundaries, we could 
consider general boundary conditions: 

f u+(0,t) = a l (t)u7(0,t)+f3 l (t), iite I ouU 
\ uT(L u t) = ai(t)u+(Li,t)+0i(t), if i G 0Ut . 

For <Xi(t) — 1 and /8»(t) = 0, we just recover the standard no-flux boundary condition 

u~l(.,t) = u^(.,t) (which is equivalent to v(.,t) = 0). 

On the outer boundaries, we also consider no-flux (Neumann) boundary conditions for <j>, which read 

d x cj> i (.,t)=0. 



(2.3) 



(2.4) 



(2.5) 

The no- flux boundary conditions mean that, on the boundary, the fluxes of cells and chemoattractants are 
null. This condition could be generalized, for example in the case when we assume that there is a production 
of fibroblasts on the boundary. 

2.2. Transmission conditions at a node 



Now, let us describe how to define the conditions at a node; this is an important point, since the behavior 
of the solution will be very different according to the conditions we choose. Moreover, let us recall that the 
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coupling between the densities on the arcs are obtained through these conditions. At node p € A/", we have to 
give values to the components such that the corresponding characteristics are going out of the node. Therefore, 
we consider the following transmission conditions at node: 

' u-(Li.t) = fijV^ii*) + E fij^i (°»*)' if * G J P' 

(2.6) 

= E t*j»t( L i>t) + E &j«7(°>*)> iH e 

where the constant £jj £ [0, 1] are the transmission coefficients: they represent the probability that a cell at a 
node decides to move from the i— th to the j — th arc of the network, also including the turnabout on the same 
arc. Let us notice that the condition differs when the arc is an incoming or an outgoing arc. Indeed, for an 
incoming (resp. outgoing) arc, the value of the function uf (resp. u~) at the node is obtained through the 
system and we need only to define u~ (resp. uf) at the boundary. 

These transmission conditions do not guarantee the continuity of the densities at node; however, we are 
interested in having the continuity of the fluxes at the node, meaning that we cannot loose nor gain any cells 
during the passage through a node. This is obtained using a condition mixing the transmission coefficients £jj 
and the velocities of the arcs connected at node p. Fixing a node and denoting the velocities of the arcs by 
Xi, i £ I p U O p , in order to have the flux conservation at node p, which is given by: 

Y,\{ut{L u t)-ur(L h t))=Y< M«i"(M)-«T(M)), (2-7) 

it is enough to impose the following conditions: 

E A -.., A,..y: !,,. (),. (2.8) 

iei P uo p 

Notice that, condition (|2.7|) . can be rewritten in the u — v variables as 

Y,v i (L i ,t)= E^(M)- (2-9) 

This condition ensures that the global mass fj,(t) of the system is conserved along the time, namely: 

fi{t) = j Ui{x, t)dx = fio := / Ui{x,Q)dx 1 for all t > 0. (2-10) 

2.3. Dissipative transmission coefficients for the hyperbolic problem. 

It is sometimes useful to restrict our attention to the case of positive transmission coefficients of dissipative 
type, in the sense that they ensure energy decay of the solutions to the linear version of system (|1.1[) . namely: 



U t + v x = 
v t + \ 2 u x = 



a.. _ ., (2-11) 



on a general network, with no-flux conditions (|2.4[) on the external nodes, and transmission conditions (|2.6I) at 
the internal nodes, always assuming the flux conservation condition (|2.8[) at nodes. 



6 



HYPERBOLIC CHEMOTAXIS ON A NETWORK 



To obtain the decay in time of the energy, which is defined by 

vl(x,t) 




E(t)= V K0M) + 



A? 




it is sufficient to impose some equalities on the coefficients, as proved in 



Proposition 1 ( [12). The energy associated with the solutions to system (|2.11l) , with no-flux conditions (I2.4[) 
on the external nodes, and transmission conditions (|2.6[) at the internal nodes, assuming condition (|2.8p . is 
decreasing if the transmission coefficients 6j belong to [0, 1], and at every node p G TV, we have: 

6j = 1 f° r al1 i£-f P U O p . (2.12) 

jei P uo p 

Actually, in |12j . it is proved that under the assumptions of Proposition [TJ it is possible to define a monotone 
generator of semigroup, and then a contraction semigroup, in the Sobolev space H 1 , for the linear transmission 
problem (|2 . 1 1 1) on a network. Let us remark also that in the simplest case of a network composed by two arcs 
(one incoming and one outgoing, see next Figure [5]) , these conditions are also necessary in order to have the 
dissipation property. In such a case we have that dissipativity is given iff: 

max jo, < 6,1 < 1, A 2 (l - 6,2) = Ai(l - 6,i)- (2.13) 

Using the previous relations and conditions on the coefficients 6j given by (|2.8[) . we obtain the values for the 
two missing coefficients: 

6,2 = 1-6,1, 6,1 = ^(1-6,0, (2.14) 
so, we have only one degree of freedom. 

2.4. Transmission conditions for 

Now let us consider the transmission conditions for <fi in system (|1.1|) . We complement conditions (|2.3j) . (|2.5[) . 
and (|2.6j) with a transmission condition for (f>. As previously, we do not impose the continuity of the density of 
chemoattractant <fi, but only the continuity of the flux at node p € J\f. Therefore, we use the Kedem-Katchalsky 
permeability condition [IS], which has been first proposed in the case of flux through a membrane. For some 
positive coefficients Kij, we impose at node 

Did n (f>i = Y K i^i ~ &)> 'e J P U °p- (2-15) 
jei p uo p 

The condition 

= Kj,i,i,j = l,..-,N (2-16) 
yields the conservation of the fluxes at node p, that is to say 

Y Did n <t>i=0. 
i£ipUO p 

Let us also notice that we can assume that Kij = 0, i = 1,...,N, which does not change condition (|2.15p . 
Finally, notice that the positivity of the transmission coefficients Kij, guarantees the energy dissipation for the 
equation for <fi in (jl.ll) . when the term in u is absent. 
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2.5. Stationary solutions 

First we consider stationary solutions, which are known to drive the asymptotic behavior of the system. Let 
us consider the case of stationary solutions of system (jl.ip . complemented with boundary conditions (|2.4[) . (|2.5p . 
(|2.6p . and (|2 . 15[) . In the general case, we find on each arc the following solution : 



Vi — constant 
u 

I hO, rr — (I, II, - b. 



exp^/Af) (Ci - -| £ eM-My)/^)dy^j , (2.17) 



which leads to solve, on each arc, the scalar non-local equation: 

-D i 4> i!XX = a i exp(<p i /\j)\C i -^ J exp(-0;(y)/Af )dy^ - bifa, (2.18) 

which has to be coupled at each node by the boundary conditions (|2.4[) , (|2.5[) , (12.61) , and (|2.15l) . 

We can prove easily that in the case of dissipative coefficients &j satisfying (|2.8[) . (|2.12p and the condition 
£ij > 0, if all the fluxes Uj are null, then the density it is continuous at a node, namely at a node p, the functions 
Ui, i S IpUOp have all the same values. However, this is not the general case. 



II 



L2 



LI 



02 



Figure 2. One incoming and one outgoing arc connected at a node. 



For the simplest network composed of one incoming / = {1} and one outgoing O — {2} arc, represented in 
Fig. [5J we find on each interval that v\ = = from condition (|2.4p . and so we obtain the following local 
system for 4>\ and 4>2 ■ 

f -Dx(j>i tXX = a\C\ exp(</>i/Af) - &i0i, f2 191 

\ -D2<p 2 ,xx = a 2 C 2 exp^/Af) - &202, 
with boundary conditions (|2.5I) and (|2. 15|) for 0i and which reads 

dMLx) = d x 4> 2 (0) = «i,2(&(0) - 0i (L0), 

and 

^i(0)=^ 2 (L 2 )=0. 
We have also to take into account the following condition given by transmission condition (|2.6p : 

A 2 6,iCi exp(0i(Li)/A?) = Ai£i j2 C 2 exp(02(O)/A!). 

Solving the corresponding system for and 02 is a difficult task, even numerically, since an infinite number 
of solutions exist both for cf>i and 2 , as in the case of a single interval [M] , and it should be necessary to make 
them verify the above conditions at node. In order to simplify our study, we limit ourselves to state a result in 
the case of constant (in space) stationary solutions to system (jl.ip . 

Proposition 2. Let us consider a general network G — (AT, A) and system (jl.ip set on each arc of the network, 
complemented with boundary and transmission conditions (|2.4p . (|2.5p . (|2.6p . and (|2.15p . 

(i) For general values of transmission coefficients satisfying (|2.8p . there is no non trivial constant 
stationary solution. 
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(ii) For the special case of transmission coefficients £ij satisfying the dissipation relations (|2.8[) and (|2.12[) 
and of the ratios ai/bi being equal to the same constant on each arc, there exists a one-parameter stationary 
solution, which is constant by arc. 

Proof. Take a constant (in space) stationary solution to system (II. ip . This means that on each arc of the 
network, we have three constant values (i^, v^, </>i), which satisfy Vi = 0, since Vi = Ui4>i X = 0, aiUi = and 
boundary conditions (|2.6[) . (|2.15p . which become in that case 

jeipUO p 

and 

= I>w(^ ~<t>i)- (2-21) 

We remark that conditions (|2.4[) and (|2.5[) are automatically satisfied. 

(i) Denoting by iV the number of arcs of the network, we have to fix therefore N unknowns to determine 
the stationary solution. Conditions (|2.20l) - (|2.21[) impose 4 equations by arc, unless the arc is connected to an 
outer node. In that case, there are only 2 conditions. To sum up, if we denote by N out the number of outer 
nodes, we need to satisfy AN — 2N out conditions. Taking into account relations (|2.8I) . we obtain that equations 
(12. 6|) are linked and the system can be reduced to a system of 4iV — 2N ou t — Ni n conditions, where Ni n is the 
number of inner nodes, which is, generally speaking, greater than the number of unknowns. Therefore, unless 
some particular sets of coefficients m,j and £i t j, the only solution for previous system is the null one on each 
arc. 

(ii) Now, let us consider transmission coefficients £i j satisfying relations (|2.8[) and (|2 . 1 2[) . We also assume 
that there exists a constant a such that, for all i, we have = abi. In that case, we can find a stationary 
solution defined on each arc by (U,0,aU). Such kind of solution satisfies clearly the transmission condition 
(|2.21j) . but satisfies also condition (|2.20l) with relations (|2.12|) . □ 

In the case (i) of the previous proposition, since the total initial mass is strictly positive and is preserved in 
time, we cannot expect the system to converge asymptotically to a stationary state which is constant on each 
arc and so non-constant asymptotic solutions are expected. In the case (ii), the constant state can be reached, 
and U is determined by the total mass of the initial data. 

3. Numerical schemes 

Here we introduce our numerical schemes. We first give some details about schemes for system (jl.ip on a 
single interval and the discretization of boundary conditions presented in [25] . Therefore, our main goal will 
be to generalize these schemes to the case of a network. In the two first subsections, we will concentrate on 
the discretization of the hyperbolic part, whereas the discretization of the parabolic part will be treated in 
subsection 13.31 

3.1. Short review of the results from |25j about AHO schemes for system (|3.1|) on a single 
interval 

Let us consider a fixed single interval [0, £]. We define a numerical grid using the following notations: h 
is the space grid size, k is the time grid size and (xj,t n ) = (jh, nk) for j — 0, . . . , M + 1, n G N are the grid 
points. In this subsection, we denote by u>™ J the discretization of function w on the grid at time t n and at 
point Xj for j = 0, ... , M + 1 and n > 0. We also use the notation / nj for f(xj,t n ), where / is an explicitly 
known function depending on (x,t). Here we describe the discretization of system (|1.1[) with no-flux boundary 
conditions v(0, t) = v(L, t) — 0, denoting by / = (f> x u and omitting the parabolic equation for (j). Since we also 
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work with Neumann boundary conditions for the <j> function, the function / will satisfy the following conditions 
on the boundary : f(0,i) — f(L,i) = 0. We therefore consider the following system 



u t + v x = 0, 
v t + \ 2 u x = f - v 

and rewrite it in a diagonal form, using the usual change of variables (|2.1[) , 

u+ + Xu+ = -(«- -« + ) + ^/- 
Set ui = ( U + ) , so that we can rewrite the system in vector form 



(3.1) 



(3.2) 



uj t + Alo x = Buj + F, (3.3) 
with A = ( ^ ? | , B = — ( ^ \) and F = — ( / | . As shown in [25] , to have a reliable scheme, 

V0Ay'2V 1 - 1 / 2A V / / 

with a correct resolution of fluxes at equilibrium, we have to deal with Asymptotically High Order schemes in 
the following form : 

, ,n+l,i , ,n,i A \ 

k 2h K ' 2h ^ 



With the following choice of the matrices 



£=-1.0.1 £=-1,0,1 

(3.4) 



o 1 / ■ 2V 01 / 2V°° 

we have a second-order AHO scheme on every stationary solutions, which is enough to balance the flux of the 
system at equilibrium. This means that the scheme is second order when evaluated on stationary solutions. 

Ah 

Monotonicity conditions are satisfied if h < 4A and k < - — ^ , see [25] for more details. Let us mention that 

it should be easy to consider third-order AHO schemes, but for simplicity (these schemes require a fourth- 
order AHO scheme for the parabolic equation with a five-points discretization for <j) x ), we prefer to limit our 
presentation to the second-order case. 

Boundary conditions for scheme Q3.4[) have to be treated carefully, to enforce mass-conservation. In [25] , the 
following boundary conditions were used : 

W n+1,0 = V n+X,M+1 = Qj 

A ft J" + V 2Ayl 2A ; ' (3.6) 

U n+1M+1 = ( X _ A\ u nM+l + ;* „,M + fc _ M y n,M + ±_ /n ,M_ 

\ h J h \h 2X J 2A 

These boundary conditions have been obtained by calculating the difference of the discrete mass at two successive 
computational times and defining u n+1 '° and U ™+ 1 > M + 1 a s a function of the discrete quantities computed at 
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time t n in order to cancel exactly this difference. Consequently, the discrete mass will be preserved in time as 
f L 

the continuous mass / u(x,t)dx is conserved for system (I3.1|) with boundary conditions v(0,t) — v(L,i) = 0, 
Jo 

at the continuous level. This technique will be generalized in this paper to the case of a network. 

3.2. The AHO scheme for system (|3.ip in the case of a network. 

Let us consider a network as previously defined in Section [21 Each arc € A, 1 < i < -/V, is parametrized as 
an interval = [0, Li] and is discretized with a space step hi and discretization points x\ for j = 0, . . . , Mj + 1. 
We still denote by k the time step, which is the same for all the arcs of the network. In this subsection, we 
denote by w™' 3 the discretization on the grid at time t n and at point x\ of a function Wi, i = 1, . . . , N on the 
i-th arc for j = 0, . . . , Mj + 1 and n > 0. 

Now, we consider the AHO scheme (|3.4[) on each interval, and we rewrite it in the u — v variables thanks to 
the change of variables (|2.ip . in order to define the discrete boundary and transmission conditions. We keep the 
possibility to use different AHO schemes on different intervals and therefore the coefficients of the scheme will 

be indexed by the number of the arc. Let R = f ^ ^ J be the matrix associated to the change of variables 
(|2.1I) , namely such that ( U ] = R I U + ] . We rewrite (13.41) in the variables u and v as : 



v \ u 



'=-1,0,1 



£=-1,0,1 £=-1,0,1 / 




(3.7) 



■7/ - 



£=-1,0,1 



£=-1,0,1 £=-1,0,1 / 

with coefficients f3 l u u , /3^ v , f3 e v u , (3 e v v and 7^, 7^ defined by 

^-'=K5c l^^-^'-K: If)- (3 - 8) 

Now, we define the numerical boundary conditions associated to this scheme. As before for equation (|3.6I) , 
we need four boundary or transmission conditions to implement this scheme on each interval. Considering an 
arc and its initial and end nodes, there are two possibilities: either they are external nodes, namely nodes 
from the outer boundaries linked to only one arc, or they are internal nodes connecting several arcs together. 
The boundary and transmission conditions will therefore depend on this feature. Below, we will impose two 
boundary conditions (|3.9[) - (|3.12[) at outer nodes, and two transmission conditions (13.10p - (|3.13[) at inner nodes. 

The first type of boundary conditions will come from condition (|2.4p at outer nodes : 

{n+1,0 n -r ■ _ T 
V, t = 0, if Z G lout, 

n+l.Mi+1 n . p • _ n 
Vi = 0, if % e Oo«t, 
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where I ou t (resp. ou t) means that the arc is incoming from (resp. outgoing to) the outer boundary. The 
second one will come from a discretization of the transmission condition (|2.6|) at node p, that is to say 



( n.Afi+l 



n,0 V* /■ ra,M,-+l V - ^ t n.O ■<• ■ _ ^ 

J 6 O p 



(3.10) 



However, these relations link all the unknowns together and they cannot be used alone. An effective way to 
compute all these quantities will be presented after equation (|3.13[) below. We still have two missing conditions 
per arc, which can be recovered by imposing the exact mass conservation between two successive computational 



N 



steps. The discrete total mass is given by 1^ ot = X^X™, where the mass corresponding to the arc % is defined 



as: 



„0 Mi 



7»,M 4 +1 



3 = 1 



Computing IgJ 1 - I™ ot , we find: 



(3.11) 



JV 



-''tot J tot 



E'HK / 1, n+1.0 n,0\ . 1 / n.l , n.0\ . -\i / n.O n,l\ . o-l n.O ol n 



}_al n,l . 1 a— 1 n,0 J_ / 1 <-n,l _ -1 <rn,o\ 





n, 












1 n 



Ln 1 / n,M 4 , n,Mi + l\ , / n,M 4 +l n,M ( \ , ol n.M; + l 



1 n,Mi 



We are going to impose boundary conditions such that the right-hand side in the previous difference is exactly 
canceled. On the outer boundaries we obtain the following type of boundary conditions, following equation 



n+1,0 



if 

r-(7«,i/j - 7»,i/i ■ ), lf * e lout, 



u 



A 

n+l,M ( +l 



1 j.n.Mi + 1 -1 j-n, Mi \ :r • _ ^ 



n,Mj , 7 | 1 Pu,v,i \ n.Mi 



(3.12) 



where J ou t and 0t ,t have the same meaning as previously. These expressions correspond to boundary conditions 
(|3.6I) in the case of a more general AHO scheme [25]. Then, using the conditions (13.121) to simplify the 



computation of X^" 1 — I^ ot and summing with respect to the nodes instead of the arcs, we can rewrite 
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the remaining difference of mass in u variables as: 



\ - \ - hik I 1 n+1 Mi + l . 1 ri+l,M s + l . n.Mi + 1/ 1 . ol . a l \ , ra,Mi + l/ 1 . „ Aj . fl l 

pf A' IS / V 

— ft 1 "\ _ ".Mi^o 1 ^ I ft— 1 I o— 1 \ I n.Mi (_o—l , o-l ^ , }_( 1 nn,M«+l _ -1 fn,Mi\ \ 
Pu,v,i) U +A \ fo, ' Pu,u,i i Pu,v,i) ' U -,i \ Pu,u,i ' Pu,v,i) > Uu,i/j ~Ju,iJi ) J' 

Therefore, using the transmission conditions (|3.10[) for U ™+ 1 ^ M «+ 1 if i £ I p and for itT^ 1,0 if z £ p , we can split 
the equation interval by interval and obtain the following numerical boundary conditions: 



n+l,M 4 +l 



-l 

A,- 



= h i [h i + E ] x (<:f - wlu,i - + u -f * +1 i 1 - 2fc f - win, 



(3.13) 



n+1,0 7 
U ■ = tli 



-1 

A 



+ K'Mm + + te -^( 2 ^ + #,«,< - + - '°)) , if i G Op- 

Once these quantities are computed, we can use equations (|3.10|) at time i n +i, to obtain u™^ 1,M * + if i £ 7 p 
and u"^ 1,0 if i £ O p . 

In conclusion, we have imposed four boundary conditions (|3.9p . (I3.10[) . (|3.12p . and (13. 13)) on each interval. 
Conditions (13.9[) and (|3.12p deal with the outer boundary and are written in the u — v variables, whereas 
conditions (13. 10)) and (|3.13p deal with the node and are written in the u ± variables. Under these conditions, 
the total numerical mass is conserved at each step. 

Now, we have to discuss the consistency of all these conditions. First, conditions (|3.9p . (13. 10)) are imposed 
exactly. Besides, it has been proved in [25 that conditions (|3.12p . set on the outer boundary, are generally of 
order one and of order two on stationary solutions. Finally, we need to consider the consistency of the conditions 
(I3.13P at node. We present here only the case i £ O p . Expanding in Taylor series up to order one, we get: 

u «+l,0 _ (1 + £ ' x ^n,0 (1 _ 2k K _ ^ _ + u „,o (l _ ^-1 . + 



3ei P uo p jei P uo p 
+ 0{k+ E hi)- 



ieipUOp 
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Now, to have consistency, namely to cancel the last two terms on the R.H.S., we need to impose the following 
condition linking the space and the time step on each arc : 

hi = 2k\i, (3.14) 

which implies, thanks to (|2.8I) : 

Under this condition and using equations (|3.2[) . expanding in Taylor series up to order three we find: 

+ few+,i(/3i )1t)i + ^ )t)ji ) + fcu" 4 l (2^- + I3l ui - Pl vt ) + Y^llJi' 1 ~ T«,i/"'°)) 
fc A; 2 

h 

_ l 2 \2 o n,0 i _ 71,0/a-l . o— 1 _ ol _ ol \ t.2 \ / /ol , ol \ q 7i,0 

K \O xx U_ i -\- ^U +i \ y p uui -\-p uvi P u ^u,i Pu,v,i) K A i\Pu,u,i Pu,v,i)°x u -\-,i 

= «,0/_ -, ol , a-1 _ al . ^ +u ™> a + fl- 1 - R 1 ■ - B 1 ) 

2 \ r^U,V,l 1 r^U,U 7 l rW,U.J "u,V,lJ 1 + V 1 t U,U.l 1 r^U^V,'l r^U.U^l r'u.V,lJ 

- h-\ + C< + - «i + + o(k 3 ). 

Thanks to this development we can state our general result of consistency. 

Proposition 3. Given a general scheme in the form (|3.4p . the conditions (|3 . 13[) at node are consistent only if 
on each arc the condition (|3.14p is verified. To have the second order accuracy at node the following conditions 
on the coefficients of the scheme have to be verified: 

@u,u,i @u,u,V Pu,v,i ~ @u,v,i ^> 7u,i — 7«,i 1- (3.15) 

Moreover, to have a third order accuracy for stationary solutions, we need : 

(3.16) 

Notice that, all these conditions are satisfied for the Roe scheme defined by (|3.5p . 
3.3. Discretization of the parabolic equation for <f> in system (|2.2p 

Now, let us explain how to compute the approximations of the function / on the arc i at discretization 

point x\ and time t n+ i needed for computing (|3.7j) . (|3.12p . and (|3.13[) . Referring to system ()2.2[) . we have 
/ = 4> x u, where <j> satisfies the parabolic equation <f) t — D <p xx = au — b<p on each arc. Boundary conditions for 
4> are given by equations (|2.5I) on the outer boundary and (|2.15p at a node. 

We solve the parabolic equation, using a finite differences scheme in space and a Crank-Nicolson method in 
time, namely an explicit-implicit method in time. 
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Therefore, we will have the following equation for f , 1 < j < Mi, 



in+l,j 



Afc 
2h? 



( + 1 I <1J. n -j J. n -j~l\ / ,71 + 1.7 + 1 . n /tl+l,j 171+1, j — l\ 



. / n+1,7 i / in-\~l,i , 



(3.17) 



Now, let us find the two boundary conditions needed on each interval. As in subsection 13.21 the boundary 
conditions will be given in the case of an outer node and in the case of an inner node. On the outer boundary, 
condition (|2.5[) for (f> is discretized using a second order approximation, which is 



,71,0 



3*< 



1 



n.2 



if i € In 



h = 3^ - 3^ , if * G O c 



(3.18) 



Let us now describe our numerical approximation for the transmission condition (|2.15p which, as the transmission 
condition for the hyperbolic part (|2.6|) . couples the <f> functions of arcs having a node in common. 

Condition (|2 . 15[) is discretized using the same second-order discretization formula as before, namely we have 
at node p, 



A n,M; + l 



■ Mi. 



■Mi-1 



1 jtip 



,n,Mi + l\ . 2 hi v 



n,0 
i,3 K^j 



,n,Mi + lN 



if i G I v 



j'eo p 



, n,0 , n. 1 

6 * = 3^ 



1 ,n.2 . 2 /lj ^ I ,n,M j+1 



2 hi yr-^ r±n,Q in,0\ - t ■ r- r\ 



3L> 



These relations can be rewritten as 



2^1 
3 Di 



jei P uO p 



.A/, 



1 ,n,M<-l 

3*< 



2 /li \ -t ,n,A/j+l 



2 r-s , n ,0 -r ■ _ r 

1 jeOp 



2_^ 

3 Di 



jeipUO P 



4 W 
3^ 



3^ + 3^ 



2 ^ V «• • f 6 n,M ^ +1 4- --^i V -^ n '° if 7 F O 



j'6-fp 



(3.19) 

Let us remark that the previous discretizations are compatible with relations (13. 18[) considering that for outer 
boundaries the coefficients Kij are null. Therefore, in this case, the value of r]° ut is just equal to 1. Since 
equations (|3.19p are coupling the unknowns of all arcs altogether, we have to solve a large system which 
contains all the equations of type (|3.17j) and also the discretizations of transmission conditions (|3.19[) . 

Once the values of J are known, we can compute a second-order discretization of the derivatives of cf> 
which gives the values of the / function, namely : 



' — /> +1 ^ +1 
2 hi V * 



b n+l,j 



1 

2h, 



7,71+1,2 



,71+1,1 



, 1 < j < Mi 



71+1,0 



),J = 0, 



1 f j.n+l.Mt—1 Aj.n+l,Mi , o ±n+l.Mi + l\ ■ jir , 1 

2^7 ^ ),j=M i + l. 
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The discretization of / needed at equations (|3.7I) . ()3.12[) . and p,13[) is therefore given by /™ + ' 3 = <^™"t ' J . 

4. Numerical tests 

Here we present some numerical experiments for system (jl.ip on networks, with the use of the methods 
introduced in Section [3J namely the second-order AHO scheme for the hyperbolic part, complemented with 
the Crank-Nicolson scheme for the parabolic part. We start with a simple test for the AHO scheme on the 
hyperbolic part of Section [3J in the case of a simplified system, where <j) x is equal to a constant a on each arc, 
for which we know the exact stationary states. 

4.1. Case <j) x constant. 

For this example, we omit the equation for <p so that the system becomes 

\u+ = -L (( a _ A)u+ + (a + \)u~) , 

1 (4-1) 
XUx = ~2A ^ a " X ^ U+ + ( a + A ) M ) • 

This system is suitable to test the accuracy of the numerical approximation, since it is easy to compute its 
asymptotic stationary solutions. We also rewrite the previous system (|4.ip using the usual change of variables 
(|2.ip which gives 

'«* + 5 = °. (4.2) 

V t + A tlx — OtU — V, 

with a a constant. To satisfy the subcharacteristic condition in [24], we also assume that 

\>\ot\. (4.3) 

Let us explain how to find the stationary states in the case of the two-arcs network of Figure [5] The method 
can be easily generalized to more complex networks. In that case, the stationary solutions satisfy the following 
equations on the intervals 1\ and li : 

Vi,x = 0, 

a}ui }X = a>i Ui - Vi, 
that is to say 

Vi = constant, 
Ui = Ci exp(aiX / Xj) + Vi/a 

Since both intervals are connected to the outer boundary, due to boundary condition (|2.4[) . we have v\ = v% = 0. 

u C ' 

Therefore we obtain non constant solutions on each arc, given by u t = = — - exp(a?i2;/A^) and the constants 

Ci are computed thanks to condition (|2.6I) . Remark that, in that case, we do not expect to have asymptotic 
states given by constant stationary solutions, since the only possible constant solution is the null one, which 
will be unsuitable, due to the constraint of the conservation of mass. Set 

Ci = ^ exp(a 1 L 1 /A?), C 2 = ^. (4.5) 

Al A2 



,„,./„,. (4-4) 



These constants solve the following system 



Al^2,l ^2(^2,2 — 1 




MC=( X ^r 1] Jl^_ u II ^ 1=0. (4.6) 
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According to (|2.8I) . Kcr M ^ {0}, and so we have at most one equation and two unknowns. Therefore, there 
exists at least one family of non trivial stationary solutions to system (|4.2I) and exactly one family when 
dim Ker M = 1. Remark that in the general case of a single node with an arbitrary number of incoming and 
outcoming arcs, assuming that all coefficients £jj are strictly positive - or more generally, that the matrix 
formed by these coefficients is irreducible, which is somewhat meaningful in the biological context -, we can 
prove that we have exactly dim Ker M = 1, thanks to the classical Perron- Frobenius theorem. 

In the case we are looking for an asymptotic state as a stationary state of the system, we can also take into 
account the conservation for mass. In that case, the stationary state we compute should have the same mass 
as the initial datum. More precisely, according to equation 

"° = tJo* Ci exp (5) dx = 5 a « ( exp ($) _ ' 

we have that the free parameter is fixed by the mass conservation. 

In particular we set L\ = 4, Li = 1, oti = a = 0.5, Ai = 2, Aa = 1 and take the dissipative transmission 
coefficients = 0.8, £2,1 = 0.4, £ 12 = 0.2, £ 2 ,2 = 0.6. If fi = 250, the system is solved by C\ ~ 28.13 and 
C2 ~ 56.25, so that the stationary solutions are u\ = C\ exp(x/8) and ui — C2 exp(x/2), with C\ ~ 34.12 and 
Ci ~ 56.25. The numerical simulations provide the asymptotic densities plotted in Fig. [3] and we notice a nice 
agreement with the stationary solutions computed analytically. Remark that densities are continuous at the 
node as explained in Section [2.51 for dissipative coefficients and vanishing fluxes. 

100 



90 
80 
70 
60 
50 
40 
30 

12 3 4 5 

FIGURE 3. Comparison between the densities of the exact and the numerical stationary 
solutions on arcs 1 and 2 obtained for Ai = 2, A2 = 1, a.% = a = 0.5, initial mass no = 250 
distributed on the network as a symmetric perturbation of the value 50, L\ = 4, L2 = 1, 
dissipative coefficients £1,1 = 0.8, £2,1 = 0.4, £i, 2 = 0.2, £2,2 = 0.6 and time T = 28. 

In Fig. |4]we present the log-log plot of the error in the L 1 -norm computed using the formula (|4.9[) of Section 
4, between the approximated and the asymptotic solutions to system (|4.2|) . The results in Fig. [4] show that the 
AHO approximation scheme provides the stationary solutions of the simplified hyperbolic model (14.21) with an 
accuracy of first order, and the error for the flux function v tends clearly to zero, faster than for the function 
u. 




HYPERBOLIC CHEMOTAXIS ON A NETWORK 



17 



10 
10" 2 

10- 3 

10" 4 

10- 5 
10- 6 
io- 7 

10" 8 
1C 

Figure 4. Log-log plot of the error in L 1 norm between the approximated and the the 
asymptotic solutions, as a function of the space step, for the solutions to system (|4.2j) . The 
error is displayed in blue for u, and in green for v. Initial data are distributed on the network 
as a symmetric perturbation of the value 50. We used different space steps satisfying condition 
(gHH , with Ai = 2, A 2 = 1, Li = 4, L 2 = 1, = 250, T = 100. 

20.6 I 1 1 1 1 1 1 1 1 

20.4 - 

S 20.2 - 

- 20 1' 

19.8 - 

19.6 I ' ' ' ' ' ' ' 

1 2 3 4 5 6 7 8 

Figure 5. Initial data corresponding to the total mass = 160. 

More examples and results showing the asymptotic behavior of solutions to the simple problem (|4.2[) on larger 
networks can be found in [2], while some analytical results are given in |12j . 

4.2. Asymptotic solutions to the full system (12.21) 

Next, we deal with the full system (|2.2p . which now include the chemotaxis equation. First, we consider 
again a network with only two arcs. We take the following data: the total mass /i = 160 distributed as a small 
perturbation of the value 20 on two arcs of length L\ — 6 and L 2 = 2, see Fig. [5J o» = 6< = 1, Ui(x 7 0) = <pi{x, 0) 
and Vi(x, 0) = 0, i = 1, 2 and Ai = 5, A 2 = 4. In the next figures we represent the asymptotic stable solutions to 
system (|2.2I) on the two-arcs network, produced by our scheme. All the solutions are plotted at a time where 
the stationary state is already reached. In particular, in Fig. |6]we plot a constant solution obtained using the 
dissipative transmission coefficients of Section 12.31 In that case we can observe what was explained in Section 
12. 5[ namely that in the case of two arcs and one node, there exist particular dissipative transmission coefficients, 
such that the asymptotic stationary solutions are constants on all the arcs. In Fig. [7] we plot the more common 
case of non-constant solutions, obtained using different parameters and non dissipative coefficients. In both 
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cases the limit flux function v is equal to zero everywhere, since for the stationary solution the flux is constant, 
the flux on the external nodes is zero, and all the arcs are connected to external nodes. 



21 
20.5 

20 
19.5 

19 

21 
20.5 

20 
19.5 

19 



~i 1 1 1 r 



_i I I I i_ 




Figure 6. Asymptotic solution for Ai 
0.25, 6,2 = 0.2, £2,2 = 0.75, T = 7.7. 



5, A2 = 4, dissipative coefficients £1, 



0.8,6 




Figure 7. Asymptotic solution at time T = 30 for Ai = 5, A2 = 4, in case of non-dissipative 
coefficients £ M = 0.8, £>a = 0.25, £1,2 = 0.24, £2,2 = 0.7. 

Let us now consider a larger network composed of twelve nodes and four arcs, see Fig. [8] We choose some 
non dissipative transmission coefficients, given in Table [TJ in order to satisfy condition (|2 .8[) . Let us consider 
as initial condition on the incoming arc 5, the function plotted in Fig. [9j where we put a small symmetric 
perturbation of the constant state u — 110. 

In this case it is hard to compute analytically the stationary solutions. We only know that non-constant 
solutions are generally expected, according to the discussion in Section [231 In Fig. [TU]we plot the asymptotic 
densities on the network node by node, starting from North-East and proceeding in a clockwise direction. Notice 
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Figure 8. A network composed of twelve arcs (six incoming and six outgoing) connected by 
four internal nodes. 





62,12 = 


= 0.1, 


61,12 = 


= 0.3, 


6,12 


= 0.3, 


6,12 


= 0.3, 


Node S-W 


62,11 = 


= 0.2, 


61,11 = 


= 0.2, 


6,11 


= 0.3, 


6,11 


= 0.3, 




62,3 = 


0.2, 


61,3 = 


0.2, 


6,3 


= 0.4, 


6,3 


= 0.2, 




6.2,4 = 


0.5, 


61,4 = 


0.1, 


6,4 


= 0.2, 


6,4 


= 0.2, 




£3,3 = 


0.1, 


60,3 = 


0.3, 


6,3 


= 0.3, 


6,3 


= 0.3, 


Node S-E 


6:10 = 


0.2, 


60,10 = 


= 0.2, 


6,10 


= 0.3, 


6,10 


= 0.3, 
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0.2, 


60,9 = 


0.2, 


6,9 


= 0.4, 


6,9 


= 0.2, 




6.2 = 


0.5, 


60,2 = 


0.1, 


6,2 


= 0.2, 


6,2 


= 0.2, 




6,1 = 


0.1, 


6,1 = 


0.3, 


6,1 


= 0.3, 


6,1 


= 0.3, 


Node N-E 


6,2 = 


0.2, 


6,2 = 


0.2, 
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= 0.3, 


6,2 


= 0.3, 
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6,8 = 


0.2, 
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= 0.4, 


6,8 


= 0.2, 




6,7 = 


0.5, 


6,7 = 


0.1, 


6,7 


= 0.2, 


6,7 


= 0.2, 




6,5 = 


0.1, 


6,5 = 


0.3, 


6,5 


= 0.3, 


6,5 


= 0.3, 


Node N-W 


6,4 = 


0.2, 


6,4 = 


0.2, 


6,4 


= 0.3, 


6,4 


= 0.3, 




6,1 = 


0.2, 


6,1 = 


0.2, 


6,1 


= 0.4, 


6,1 


= 0.2, 




6,6 = 


0.5, 


6,6 = 


0.1, 


6,6 


= 0.2, 


6,6 


= 0.2. 



Table 1 . Transmission coefficients used for the numerical simulations of Figures [10] and [TTJ 
given node by node. 



that most of the arcs are repeated in the different figures. In Fig. [TTJ the asymptotic fluxes are represented, 
and again our scheme is able to stabilize them correctly. We notice that the fluxes of arcs connected to outer 
boundaries vanish, whereas the fluxes of inner arcs, even if they are constant, are different from zero. 

4.3. Instabilities: the appearance of numerical blow-up 

Let us consider some cases that present a strong asymptotical instability. Indeed, for some values of the 
parameters of the problem, namely of the arc's length L and the cell velocity A, in connection with the total 
mass distributed on the arcs of the network, we can observe increasing oscillations, which eventually may cause 
the blow-up of solutions. It is important to notice that the blow-up can be already observed for this model 
even for a single arc, see Example 1 below, when the total mass no is large with respect to the characteristic 
parameters L and A. However, here the presence of more arcs, and so, a greater total length and total mass, 
makes this kind of phenomenon much more frequent. 
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Figure 9. Initial condition for u and <j> on arc 5 of the network presented in Fig. [8] 




FIGURE 10. Stationary solutions for the network composed of 12 arcs and 4 nodes of Fig. 
[8] the densities are computed at time T = 30, the values of the parameters are given by: 
A, = A = 10, Li = 1, Ot = hi = Di = 1. The transmission coefficients can be found in Tabic Q] 
The total initial mass no = 1320 is distributed as a perturbation of the constant state 110 on 
arc 5 as in Fig. [9] and as the constant density 110 on the other arcs, with hi = h = 0.01, k = 
0.0005. 



Example 1. Here we assume that we have only one interval with L = 1 and A = 10 and we take, as initial 
condition for the density and the chemoattractant, a symmetric perturbation of a constant state Co = 9000. 
The total mass is no = 9000, as shown in Figure [HJ The solution presents a clear blow-up at time T = 0.1, 
see Fig. [131 This blow-up seems associated to non physical negative values of the density function u, and it is 
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Figure 11. The asymptotic fluxes of the arcs of the network composed of 12 arcs and 4 nodes 
at time T = 30, with the same data as Fig. [TUJ 



observed in the same way even for refined meshes (see Table 2 for the case of two arcs). This is not surprising, 
since the quasimonotonicity of the system, see again |24j . is violated when the gradient <f> x is larger than A. 




Figure 12. The initial condition uq(x) is a symmetric perturbation of a constant state 
Co = 9000, the total mass is fi = 9000. 



Example 2. Here we take two arcs of length Ly = 6 and L2 = 2 and the initial density as in Fig. [5j with 
aj = bi = 1, Ui(x,0) = (f)i(x,Q), and Vi(x,0) = 0, i = 1,2. Then we change the values of velocities Ai and 
A2 in order to see how they influence the behavior of solutions to system (jl.ip . At the junction we assume 
transmission and dissipative coefficients, taking £1^ ~ 0.96 and then satisfying equations (j2.13[) - (|2.14[) . What 
we observe is that solutions blow up in finite time or not according to the relative values of Ai and A2, as it is 
shown in Figure [T4l More precisely, we can observe three different regimes. If A2 is large with respect to A 1 _ 2 , 
solutions stay bounded and converge to stationary solutions (green "x" in Figure fl4|) . If Ai is small with A2 
large enough, then solutions blow up in finite time (red " +" in Figure I14p . Finally, there is a small region in 
between, Ai around the value 3 and A2 small enough, such that solutions present a large spike at the boundaries 
(marked by blue asterisks "*"). 
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Figure 13. Blow-up of the solution at time T = 0.1, for data in Fig. [12] with L = 1, A = 10, 
h = 0.001, /io = 9000: on the left the blow-up density u and the concentration (j>, on the right 
the flux v. 
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Figure 14. Schematization of the regions describing the behavior of solution for fio =160 and 
the velocities Ai and A2 varying: blow-up (marked by red crosses "+"), solutions with a spike 
at the boundaries (marked by blue asterisks "*") and stable stationary solutions (marked by 
green "x"). 



Let us now focus on the blow-up behavior. Referring to Fig. 1141 we can choose a pair of velocities belonging 
to the blow-up region marked by red crosses "+", to say Ai = 1 and A2 = 2. The time step just before the 
numerical blow-up time of corresponding solutions, starting from initial data as in Fig. is plotted in Fig. OH 
Even if apparently we are close to the transmission point, there are many grid points separating it from the 
blow-up point. To show that the blow-up is not just a numerical artifact, we perform the same simulation with 
the same data, but on refined grids. In Table [2] we report the blow-up time of solutions to system (jl.ip for a 
fixed global mass /j,q when either the CFL condition v = £ A or h go to zero. Out of the case of v = 1, which 
appears to be more unstable, the blow-up time is independent of the meshes and has to be considered to occur 
in the analytical solutions. 
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FIGURE 15. Blow-up at time T = 4, for initial data as in Fig. with L\ = 6,1,2 = 2, Ai = 1 
and A2 = 2, dissipative coefficients with £x,i = 0.96, the total mass is equal to /iq = 160: on the 
left the density u and the concentration 4>, on the right the flux v. The space steps are equal 
to hi = 0.001, h 2 = 0.002. 







Blow-up 


time 




h 


v = 1 








0.01 


2 


4 


4 


4 


0.0025 


1 


4 


4 


4 


0.001 


0.5 


4 


4 


4 



Table 2. Blow-up times of the solutions to system (|l.l|l when either the CFL condition v = \\ 
or h go to zero, with transmission coefficients of dissipative type, L\ = 6, = 2, Ai = 1, A2 = 2, 
Ho = 160. 



4.4. Comparisons and errors 

Let us now introduce the formal order of convergence of a numerical method j w for the computation of the 
solution w as the minimum among the orders on the arcs of the network: 



7,„ = nun 7,, 



(4.7) 



where 



7w = lo §2 



»=!,. 



(4.8) 



The L -error for the numerical solution on each arc is 



hi 
n 



E 



l=0,...,nMi 



W, 



n = l,2, 



(4.9) 
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where wj(h) denotes the numerical solution obtained with the space step discretization equal to h, computed 
in xj at the final time T. The total L 1 -error is 

N 

TOT^^^e'^). (4.10) 

i=l 

Table [3] shows the L 1 -error (14.91) on the asymptotic solutions u, <j> and v and order of convergence (14.7[) of 
the approximation scheme applied to the considered network. 



h 




Error on u 


70 


Error on <f> 




Error on v 


0.025 


0.916393 


1.78849e-04 


0.965238 


1.78848e-04 


1.212334 


3.34559e-07 


0.0125 


0.959614 


8.87206e-05 


0.982631 


8.87207e-05 


-0.058657 


1.44060e-07 


0.00625 


0.980243 


4.41941e-05 


0.990856 


4.41954e-05 


0.666605 


1.49949e-07 


0.003125 


0.986317 


2.20550e-05 


0.992983 


2.20651e-05 


0.863690 


9.43741e-08 


0.0015625 


0.937936 


1.10172e-05 


0.937109 


1.10280e-05 


0.955806 


5.17981e-08 



Table 3. Orders and errors of the approximation scheme for the solutions to system (|1,1[) . 
Li = 1, Ai = 4, i = 1, 2, no = 120.056, T = 25. 

The results in Table[3]show the effectiveness of AHO approximation scheme in the solution of the transmission 
problem represented by the hyperbolic model (jl.ip . We notice indeed that even in this more general case the 
scheme still keeps a formal accuracy of first order, although the interactions at the boundaries could deteriorate 
its accuracy. 
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